Structural and dynamic properties of eutectic mixtures based on menthol and fatty acids derived from coconut oil: a MD simulation study

The structural and dynamical properties of the binary mixture of Menthol (MEN) and Fatty acids (FAs) were investigated using molecular dynamics simulations. To this end, the relationship between the structural and dynamical properties of the eutectic mixtures of MEN and FAs with different molar percentages of FAs are studied. Structural properties of the eutectic mixtures were characterized by calculating the combined distribution functions (CDFs), radial distribution functions (RDFs), angular distribution functions (ADFs), hydrogen bonding networks, and spatial distribution functions (SDF). Additionally, our Results indicated robust interactions between menthol and Caprylic acid molecules Finally, the transport properties of the mixtures were investigated using the mean square displacement (MSD) of the centers of mass of the species, self-diffusion coefficients and vector reorientation dynamics (VRD) of bonds. Overall, our simulation results indicated that intermolecular interactions have a significant effect on the dynamic properties of species.


Computational methods
Three acids such as Caprylic acid (CAP), Decanoic acid (DEC), and Lauric acid (LUA) were selected in the saturated fatty acids for MD simulation. PACKMOL package was used for the preparation of the binary mixtures of MEN and FAs with the different molar percentages of FAs 11 . The information about the simulated boxes and their different combinations (HBD and HBA) with abbreviations are listed in Table 1. The initial simulation box included with the mole percent of FAs at 20.0% and 80.0% of MEN was provided using random distribution. In the next step, the number of FAS molecules was increased for the preparation of the mixtures with the mole percent of FAs at 70.0%. The NAMD 2.13 software package was used to perform all the MD simulations 15 . First, all the initial configurations were minimized to remove the bad contacts and then the binary mixtures were heated to 323 K to prepare a uniform liquid phase. In the next steps, an NVT ensemble (constant number of particles, volume, and temperature) was performed for 5 ns at 323 K which is followed by the equilibration of the binary mixture for 50 ns in the NPT ensemble (constant number of particles, pressure, and temperature) with a timestep of 1 fs. After equilibrating the systems, its liquid phase was analyzed at 353 K. The binary mixtures with the most intermolecular interaction between HBD and HBA were selected, and their temperature was increased slowly from 0 to 353 K and then decreased back to 323 K and 300 K with a temperature gradient of 1 K/FS. The Nose-Hoover thermostat 16 and the Parrinello-Rahman pressure coupling 17 were used to maintain a constant www.nature.com/scientificreports/ temperature and pressure, respectively, during the simulation. Periodic boundary conditions were imposed or the simulated systems. Force field parameters and partial atomic charges were calculated by fitting the RESP computed at the MP2/6-31G* level 18 . The applied total potential energy has the following different terms: This force field includes bonded parts such as bonds, angles, and torsional dihedral interactions and the non-bonded interactions consisting of van der Waals (vdW) and electrostatic (Eel) terms 18 . The particles mesh Ewald (PME) method and the Lennard-Jones (LJ) 6-12 potential with a cut-off distance of 12 Å were applied for long-range and short-range interactions, respectively 18 .

Results and discussion
Radial distribution functions (RDFs). The radial distribution function (RDF) can represent how the average species density varies as a function of distance from a given species within the system 18 . This quantity was computed to determine the average distribution of species around any given species the 18 . g(r) is given by where N indicates the number of species in the system, ρ is the density of species 18 . The angle brackets are an average over all the species. This information could be used to characterize various interactions between species in the binary mixtures of MEN and FAs. To better understand the structural properties of the binary mixtures, RDFs of the H atoms of HBD around the O atoms of the HBA was obtained. Figure 1 depicts labels of the atoms in the hydrogen bond donor and acceptor. The Fas-MEN RDFs have a main peak and a small plateau in the vicinity of the first maximum peak, while all RDFs between menthol molecules have a similar shape with two maxima. The sharp peaks of HA Fas -OA MEN RDFs are clear indications of interaction between the hydroxyl hydrogen of FAs and the O atom of MEN at 2 Å distance (see Fig. 2a). The pair interaction of two species has more intensity in the binary mixtures of MEN and with % CAP = 44, % DEC = 35, and % LUA = 25 (see Fig. 2b,c). The RDFs between menthol molecules were investigated in the binary mixtures with different percentages of FAs molecules. RDFs between menthol molecules in the binary mixtures with % DEC = 20 and 35 are shown in Fig. 2d. The shoulder and the main peak are located at about 3 and 2 Å in the binary mixtures of DEC and MEN. The height of the first peak of the MEN-MEN RDF gradually was decreased with increasing the mole percent of FAs. These results indicate that the probability distribution of menthol molecules around FAs molecules increases in the binary mixtures.    Fig. 4. In general, the CNs are reduced with increasing temperature. The coordination number of CT12 around menthol was notably greater than the CT4 in the binary mixtures of MLA. This aggregation is the result of the van der Waals force and it can be also due to lack of many hydrogen Coordination number where σ, and X are the standard deviation and the average number of hydrogen bonds, respectively 22 . The average number of hydrogen bonds between species is an important parameter for the formation of eutectic solvent that changes under the influence of the relative percentage of species and temperature. The hydrogen bonds in the binary mixtures of MCA, MDA, and MLA with %FAs = 44, 35, and 25 are more than hydrogen bonds in the mixtures with other ratios (see Supplementary Table S2). The H-bond number between FAs and MEN was reduced in the binary mixtures with a higher percentage of FAs. Most likely, Decreasing the average number of hydrogen bonds is due to the distribution of the number of hydrogen bonds between the carboxyl groups of FAs (see Fig. 5a). The hydrogen-bonding network is directly affected by the temperature of the mixture. Investigation of the hydrogen bonding network in the binary mixture at the temperature range of 300 to 353 also showed that the H-bond number was slightly decreased with raising of the temperature (see Fig. 5b). The total time that a hydrogen-bond network has established between HBD and HBA is described using the percent occupancy 23 . This analysis has depicted the stability and persistence of hydrogen bonds between the possible HBD and HBA pairs 24 . The stability of the H-bonds interactions between different atoms of the MEN molecules and the O atom of FAs molecules was investigated in the binary mixture. The hydrogen-bonded network between the HA atom of MEN and the O atom of FAs is more stable than the other atoms of MEN. Furthermore, interactions between the HA atoms of FAs and the OA atom MEN have the highest percent occupancy of hydrogen bonding. The stability of the hydrogen bonding network between the different species is provided in Fig. 6a,b. Figure 6a displays the Hydrogen bond occupancy percentage of the HA atoms of MEN and the O atom FAs in the binary mixture with a molar percentage of 70% at 323 K. The percent occupancy of the hydrogen bonding between HBD and HBA was the highest in the binary mixture of MCA. The stability of the H-bond between the two species was decreased by increasing the carbon chain length of FAs from C8 to C12. The occupancy of the hydrogen bonding between three kinds of HBA and MEN has the highest content in the binary mixtures of MCA, MDA, and MLA with %FAs = 44, 35, and 25, respectively, the number of H-bonds between different species (HBD and HBA) were measured over time. We surmise that, the strong interaction between FAs and MEN at 300 K is due to hydrogen bonding with the highest percentage of occupancies. The stability of the hydrogen bonding network confirmed density reduction trend of the mixture at distinct temperatures (see Supplementary Table S4).    www.nature.com/scientificreports/ distribution around the FAs molecules in the binary mixtures is more likely associated with the accumulation of acid molecules around each other. The maximum distribution between the two species of MEN and FAs is not generally the same at different temperatures. Increasing temperatures cause a shift in the aggregation behavior of molecules around each other. At high temperatures, it appears that the stability of intermolecular hydrogen bonds decreases. The stability of intermolecular hydrogen bonds may affect the local solvation structure in the binary mixtures. It should be noted that the isovalues were selected to reflect the completion of the first solvation shell 26 .

Tracing [MEN] [FAs] H-bonds using combined distribution function (CDF). The distance criteria
of the H-bond were determined from the first minimum of the peaks of the corresponding RDFs in the binary mixtures. However, the distance criteria may not be adequate to capture the presence of a strong H-bond in the binary mixtures 21 . Therefore, the combined radial/angular distribution functions were obtained for all possible pairs. The composed CDFs of the RDF between the hydroxyl hydrogen (HA) site of MEN and the O atom of the FAs as the X-axis, the Angular Distribution Function (ADF) between two vectors (R1 and R2) as the Y-axis are shown in Fig. 8a. H-bonds (weak or strong) were found at the angles range between 130° and 150° and distances of less than 2 Å. Figure 8a shows that there are more menthols around CAP than the other acids in the binary mixtures (see Fig. 8b

Non-bonded interaction energy. The non-bonded interaction energy is an important analysis for
understanding intermolecular interactions in the binary mixtures. The non-bonded energy includes two terms, short-range van der Waals (vdW) interactions and long-range electrostatic interactions, which are calculated using Lennard-Jones (12-6) function and Coulombic equation, respectively 27 . The vdW and Coul interactions are described by the following equations: The values of the non-bonded energies are reported in Supplementary Table S3. It could be seen from Supplementary Table S3, Ecoul and E vdW are − 3235.06 kJ/mol and − 8.10 kJ/mol, respectively, in the binary mixtures CAP and MEN with %FAs = 44 at 323. These results point out that the electrostatic interaction energy is predominant between HBA and HBD. The vdW energy values are so trivial that they can be disregarded. The absolute value of the non-bonded energies shows that the interaction between LUA and MEN was significantly greater in the binary mixture with %FAs = 25 at 323 K than the other ratios. In term of energy, this percentage combination may be more suitable for the preparation of the eutectic mixture based on acid and menthol. However, interactions between two molecules reveal that temperature can be a more effective factor to be considered. Investigation of the interaction energy of MEN: MEN in the binary mixture with %FAs = 44 showed that the amount of energy was slightly decreased under various temperatures from 300 to 353 K (see Fig. 10). According to Supplementary Fig. S1, increasing the length of the alkyl chain of FAs may prevent further interaction between the HBA and HBD. The reduction of non-bonded energy is most likely related to the reduction of the partial charge on oxygen atoms FAs.
Dynamical and transport properties. The mean-squared displacements (MSDs) of species were measured to describe the migration of species 28 . The MSDs of species as a function of time were plotted for binary mixtures of FAs and MEN in the same ratio, see Fig. 11a,b. The low slope of the MSD curve of menthol molecules indicates a strong interaction between the two components in the binary mixture of menthol and Caprylic acid compared to the mixture of MDA. In comparison to MSD of DES molecules, due to interactions between FAs molecules which FAs molecules move in binary mixtures, the slope of the MSD curve of LUA molecules shows smaller value. The self-diffusivity coefficient is the macroscopic dynamical property which is calculated from MD simulation 29 . The self-diffusion coefficients of species were obtained from the slopes of the lines fitted to the MSD curves using the Einstein relation: www.nature.com/scientificreports/ where ��|r(t)| 2 � is the mean-square displacement (MSD) of species i. To ensure the measured values of selfdiffusion coefficients, the diffusive regime was identified and sampling was performed in this regime 30 . The β parameter was used to determine the location of the diffusive regime in the binary mixtures 31 . The beta parameter was used to the self-diffusivities of the species in DESs by Colina and coworkers 18 and it is calculated according to: The self-diffusion coefficient of the species was derived from the slope of the log-log plots of MSD in the range of 40 to 50 ns 18 . The calculated self-diffusion coefficient values from the diffusive regime are shown in Table 2. The self-diffusion coefficient obtained for the CAP was 0.92 Å 2 ns −1 and 10.29 Å 2 ns −1 , respectively, in the binary mixture with %FAs = 44 at 300 K and 353 K. In general, the reduced viscosity at higher temperatures has led to more rapid migration of species in the diffusive regime. The self-diffusion coefficients of species in the binary mixtures with a molar percentage of 70% are ranked as MDA > MCA > MLA. The self-diffusion coefficient values for the CAP molecules were 4.29 Å 2 ns −1 and 3.25 Å 2 ns −1 , respectively, in the binary mixture of FAs and MEN with %FAs = 20 K and 44 K, at 323 K. The self-diffusions of species decrease by adding the acid molecules into the simulation box. It should be noted that the remarkable changes were found in the binary mixtures of MCA, MDA, and MLA with %FAs = 44, 35, and 25. Reducing self-diffusion can be justified by intermolecular interactions.

Vector reorientation dynamics (VRD). The vector reorientation dynamics is a common quantity that
can be determined by both experimental techniques 32 and molecular dynamics simulation methods. VRD is computed from Eq. (9), where − → a i (t) and − → a i (t + τ ) are a vector selected on a molecule in step t and at the later time t + τ, respectively 33 . The reorientation dynamics of the bonds vector of the MEN molecules is displayed in Fig. 12a. As one can see, the reorientation of the bond vectors OA-HA in MEN is slower than the other bond vectors due to intermolecular interactions. As a result, the reorientation dynamics of HA-OA vectors is discussed in this section. For the reorientation of the HA-OA bond vectors, the blue and red solid lines and the green dashed lines indicate VRD in the binary mixtures of MCA, MDA, and MLA, respectively (see Fig. 12b). VRD of the HA-OA bond in Capric acid is considerably faster than this bond in the other acids. VRD of the HA-OA bond of Capric acid in the binary mixture of MCA with %FAs = 44 was the lowest one and increased at 353 K (see Fig. 12c). Interactions between FAs and MEN resulted in the lowest orientation velocity in the mixture of MCA with %CAP = 44. Similarly, the HA-OA bond has the fastest reorientation dynamic in the binary mixture with %LUA = 25 at 353 K compared to other temperatures. The slower orientation of the HA-OA bond can be related to the interaction of the hydrogen bond between the acid molecules that are reported in the "Hydrogen bond analysis". Here, p x i and u ij indicate the momentum component and the interaction potential, respectively, and x i is a component of the radius vector 36 . The measured viscosity values are collected in Table 3. At 323 K, the values of η of the MCA and MDA binary mixtures are 121.6 and 104.51 MPa, respectively. Based on the results of the      www.nature.com/scientificreports/ occupancy of the hydrogen bonding, the hydrogen bond network in the binary mixture of CAP and MEN is more stable than in the binary mixture of MDA. Thus, it can be concluded that the interactions between FAs and MEN molecules have an important effect on the shear viscosity of the binary mixtures. The effect of increasing FAs molecules on the shear viscosity of the binary mixtures was investigated by measuring the η at the different molar percent of FAs and the results reveal that adding of the FAs molecules into the binary mixtures increases the viscosity of all the mixtures of FAs and MEN due to bigger value of η in the binary mixtures containing 44% CAP is higher (52.451 MPa) than the binary mixture with the molar percentage of 20% (24.63 MPa). The values of the shear viscosity are in reasonable agreement with the hydrogen bonds result; therefore, the shear viscosity of the binary mixtures was decreased gradually during increasing of the temperature. Finally, It should be noted that the shear viscosities obtained for the binary mixtures of LUA and MEN with %FAs = 25 using MD simulation are close to the experimental value, with a difference of less than 2% 8 .

Conclusion
In the present work, MD simulations was exploited to investigate the correlation between the structural and dynamical properties of hydrophobic eutectic mixtures. There is a strong interaction between FAs and MEN molecules in the binary mixtures of MCA, MDA, and MLA with %FAs = 44, 35, and 25, respectively. The maximum distribution between the two species of MEN and FAs in these ratios is estimated from structural analyzes such as RDF, CDF, SDF, and the density distribution function. The distance range of 2 to 3 Å and the range of angles between 130° and 180° of CDFs shows the maximum number of hydrogens bonds between the two species. A single fatty acid molecule can participate as a bridge in the binary mixture because it can accept H-bonds using the -COOH group. The increased temperatures lead to an increase in the migration of molecules, which we surmise that it is attributed to the changes in the intermolecular interactions. Furthermore, it was observed that the interaction between two species was reduced with the increase in the length of alkyl chain.